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Abstract 

We have simulated, showcase, the pentapeptide Met- 

enkephalin (Tyr-Gly-Gly-Phe-Met) to visualize the energy landscape 
and investigate the conformational coverage by the multicanonical 
method. We have obtained a three-dimensional topographic picture of 
the whole energy landscape by plotting the histogram with respect to 
energy(temperature) and the order parameter, which gives the degree 
of resemblance of any created conformation with the global energy 
minimum (GEM). 

Keywords: Energy Lanscape, Conformational Sampling, Multicanon- 
ical Simulation. 



1 Introduction 

Biological macromolecules such as proteins have a well defined 3D structure 
which is essential for their biological activity. Therefore, predicting the pro- 
tein's structure by theoretical/computational methods is an important goal 
in structural biology. |L| 

The configuration space of peptide's and protein's presents a complex energy 
profile consisting of tremendous number of local minima; their basins of 
attraction were called localized microstates. The energy profile also contains 
larger potential energy wells defined over wide microstates (e.g., the protein's 
fluctuations around its averaged structure), each including many localized 
ones. f||] 

Because of energy barriers, the commonly used thermodynamic simulation 
techniques, such as the Metropolis Monte Carlo (MC) || and molecular dy- 
namics (MD) are not very efficient in sampling a rugged landscape. Thus, 
the molecule remains in its starting wide microstate or move to a neighbor 
wide microstate, but in practice will hardly reach the most stable one. The 
system may occur to be trapped in a basin for a long time, which results in 
non-ergodic behavior. Therefore, developing simulation methods that lead to 
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an efficient crossing of the energy barriers has been a long standing challenge. 

The topography of the energy landscape, especially near the global mini- 
mum is of particular importance, because the potential energy surface defines 
the behavior of the system. Methods for searching energy surfaces are pro- 
posed ||, energy landscape perspectives are investigated || and the fractal 
dimensions are studied 0. The essence of a funnel structure of energy land- 
scape at some fixed temperatures has recently been shown by Hansmann 
and Onuchic ||. Consequently, a visualization of the whole rugged land- 
scape covering the entire energy and temperature ranges would be helpful to 
develop methods allowing one to survey the distribution of structures in con- 
formational spaces. Such a goal can be achieved within the multicanonical 
ensemble approach. 

An ideal simulation scheme should freely visit the entire configuration space 
and predominantly sample the significant conformations. The trapping prob- 
lem of the MC and MD methods can be alleviated to a large extent, by the 
multicanonical MC method (MUCA) 10, |TTJ , which was applied initially 



to lattice spin models and its relevance for complex systems was first no- 
ticed in Ref. | IU|| . Application of the multicanonical approach to peptides 
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was pioneered by Hansmann and Okamoto \T2\ and followed by others [I3"|; 
simulations of protein folding with MUCA and related generalized ensemble 
methods are reviewed in Refs. [1T4| and [f5f . 



2 The Model 

The multicanonical ensemble based on a probability function in which the 
different energies are equally probable. However, implementation of MUCA 
is not straightforward because the density of states n(E) is unknown a priori. 
In practice, one only needs to know the weights u, 

w(E) ~ l/n(E) = exp[(£ - F T{E) )/k B T{E)). (1) 

These weights are calculated in the first stage of simulation process by an 
iterative procedure in which the temperatures T(E) are built recursively to- 
gether with the microcanonical free energies Ft(e)/^bT(E), up to an additive 
constant. The iterative procedure is followed by a long production run based 
on the fixed w's where equilibrium configurations are sampled. Re- weighting 
techniques (see Ferrenberg and Swendsen and literature given in their 
second reference) enable one to obtain Boltzmann averages of various ther- 
modynamic properties over a large range of temperatures. 
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As pointed out above, calculation of the a priori unknown MUCA weights is 
not trivial, requiring an experienced human intervention. For lattice models, 
this problem was addressed in a sketchy way by Berg and Celik |TO and 



later by Berg [ l7| . An alternative way is to establish an automatic process 
by incorporating the statistical errors within the recursion procedure. The 
automatic procedure was tested successfully ]TJ] as applied to models of the 
pentapeptide Leu-enkephalin (H-Tyr-Gly-Gly-Phe-Leu-OH) described by the 



ECEPP/2 potential energy function [IS 



In this work, as in our previous one, [18| Met-enkephalin is modeled by 
the ECEPP/2 potential, which assumes a rigid geometry, and is based on 
non-bonded, Lennard- Jones, torsional, hydrogen-bond, and electrostatic po- 
tential terms with the dielectric constant e = 2. This potential energy is 



implemented into the software package SMMP ||21| . We further fix peptide 
bond angles uj to their common value 180°, which leaves us with 19 dihe- 
dral angles as independent degrees of freedom ( rip — 19 ). We have also 
simulated Met-enkephalin with variable peptide bond angles, for which the 
distribution of conformations are included in Table I. 
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3 Results and Discussions 

We first carried out canonical (i.e., constant T) MC simulations at the rel- 
atively high temperatures and MUCA test runs which enabled us to deter- 
mine the required energy ranges. Then we preformed full simulation which 
cover the high temperature region up to T max = 1000 K reliably. The en- 
ergy range was divided into 31 bins of lkcal/mol each, covering the range 
[20, —11] kcal/mol. The lowest energy encountered was —10.75 kcal/mol and 
T max = 1000 K was also used above 20 kcal/mol. At each update step, a trial 
conformation was obtained by changing one dihedral angle at random within 
the range [—180°; 180°], followed by the Metropolis test and an update of the 
suitable histogram. The dihedral angles were always visited in a predefined 
(sequential) order, going from Tyr to Met; a cycle of iV MC steps (iV = 19) 
is called a sweep. The weights were built after m = 100 recursions during 
a long single simulation, where the parameters b{ and were iterated every 
5000 sweeps. 

For peptides it is not only of interest to obtain thermodynamic averages 
and fluctuations at different temperatures but also to find the most stable 
regions in conformational space populated by the molecule. In the organic 
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chemistry community conformational search methods have been developed 
and attempts have been made to find the global energy minimum (GEM) 
and all the energy minimized conformations in certain energy ranges above 
the GEM (see Ref.2(b), and references cited therein). 

The lowest energy conformation (our suspected GEM) was found at E — 
-10.75 kcal/mol. 

Here we define, following Hansmann et. al. ||, an order parameter (OP) 

o^-5^glM s, l. P) 

where a\ RS ^ ve a® are the dihedral angles of the reference state (which 
is taken as GEM) and of the considered configuration, respectively. The 
difference af — <y\ RS ^ is always in the interval [—180°, 180°], which in turn 
gives for peptides 

< < OP > T < 1 (3) 

Figure. 1 shows the energy landscape obtained by the multicanonical simula- 
tion run of one million sweep plotted against energy and the order parameter. 
Here, we would like to point out that the utilized data is obtained by sam- 
pling of the conformational space and no minimization procedure is applied. 
At high temperatures, where the peptide is in the random coil state, the his- 
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togram looks as one gaussian-like peak centered around the value of the order 
parameter OP ~ 0.3. When the temperature is lowered, first a transition 
from the state of random coil to globular structure is expected. In Figure. 2 
we show the same energy landscape of Fig. 1(b) by grouping the conforma- 
tions of 1 kcal/mol interval in energy. Curve a) denotes the energy inter- 
val — 1 kcal/mol < E < kcal/mol, which corresponds after re- weighting 
to the temperature interval 315 K < T a < 330 K. At this temperature, 
the energy landscape starts deviating from a smooth surface and develops a 
shoulder. We identify this temperature as the starting of forming a struc- 
ture rather than a random coil. Further down in energy (temperature), the 
newly developing branch becomes more populated. At the temperature 215 
K < T& < 230 K denoted by the curve b), the energy landscape displays a 
typical structure bifurcating into two branches of almost equal height. From 
there on, the branch having larger values of the order parameter wins and 
more conformations populate that section of the conformational space. Our 
estimate of T a and Tb from the topographic structure of the potential energy 
surface of Met-enkephalin are very close to the values of the collapse tem- 
perature Tg = 295 ± 20 K and the folding temperature Tf = 230 ± 30 K, 
respectively, determined by Hansmann et al |2(| . We observe a third temper- 
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ature denoted by the curve c) in Fig. 2 where the glassy behavior sets in and 
many valley structure of the landscape become clearly pronounced. For our 
simulated peptide sample Met-enkephalin, this temperature is in the range 
155 K < T c < 185 K. Below this temperature, one observes the appearances 
of multiple valleys which are well separated by high energy barriers. The val- 
ley at the far-out end of the order parameter scale having the conformations 
with the value of the order parameter in the range 0.98 < OP < 1 contains 
the global energy minimum (GEM), respect to which the order parameter is 
evaluated. The temperature T c seems to correspond to the glass transition 
temperature estimate of T g = 180 ± 30 K, which value is based on the frac- 
tal dimension estimates. |7J In Fig. 3 we plotted all the conformations found 
with energy E < —10.5 kcal/mol with respect to the order parameter. Their 
number is 3587 conformations in one production run of one million sweeps. 
As clearly seen from Fig. 3 that the conformations in this energy range are 
localized in one of the four valleys, which are identified by the value of their 
order parameter OP ~ 0.80, 0.87, 0.92 and 0.98. The conformations in the 
neighborhood of the GEM take place within the same wide microstate of 
the GEM but they are grouped into local microstates, each of which are 
one of the above mentioned valleys. The small differences in values of OP 



8 



comes from the differences in side-chain angles. We observe no conformation 
anywhere outside the definite valleys when the energy is less than about f 
kcal/mol above the GEM. 

The number of conformations found in energy bins of 1 kcal/mol, which were 
plotted in Fig.2, appear in Table I. The lowest bin is 0.75 kcal/mol and 
includes the GEM. The table displays the distribution of sampled conforma- 
tions according to the order parameter values, namely the distribution with 
respect to how far they are in configuration space from the global energy 
minimum. We also included in Table I the same distribution obtained in our 
simulation of Met-enkephalin for the case of variable peptide-bond angles u>. 

In conclusion, we have simulated the pentapeptide Met-enkephalin by utiliz- 
ing the multicanonical ensemble approach and investigated the structure of 
the rugged energy landscape in the configurational space. We were able to 
display the distribution of at all temperatures from a single simulation and 
estimate the critical temperatures such as the collapse temperature, the fold- 
ing temperature and the glass transition temperature. Such a visualization 
would be helpful in designing algorithms for efficient sampling of conforma- 
tional space. 
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Table 1: Number of conformations in energy bins of 1 kcal/mol. 



ENERGY 


OVERLAP 


TOTAL 


Fix uj 


1.0-0.9 


0.9-0.8 


0.8-0.7 


0.7-0.6 


CONF. 


i n 7x. +n inn 
-1U. /o to -1U.U 


3282 


3935 


3073 


2779 


loo/ 1 


1 n n +n q n 
-iu.u to -y.u 


1001 


3530 


4925 


4475 


osnsg 
zoUoo 


q n fn s n 
-y.u to -o.u 


467 


2332 


4003 


3979 


zozzu 


q n tn 7 n 

-O.U LO - 1 .U 


190 


1460 


3150 


3488 




-7.0 to -6.0 


90 


897 


2515 


3290 


22497 


Variable u 






-12.21 to -12.0 


23 25 - - 


48 


-12.0 to -11.0 


6380 


7568 


302 


197 


14457 


-11.0 to -10.0 


7600 


21199 


4775 


2784 


37107 


-10.0 to -9.0 


2700 


9956 


3959 


3456 


28430 


-9.0 to -8.0 


600 


3107 


2390 


3137 


2644 
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Figure 1: Energy surface in configuration space of Met-enkephalin viewed 
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from different angles. 
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Figure 2: Same as Fig. 1(b), plotted by grouping the conformations of 1 
kcal/mol interval in energy. 
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Figure 3: Distribution of microstates with E < —10.5 kcal/mol with respect 
to the overlap parameter. 



17 



